ORIENTATIONAL PHASE TRANSITION IN SOLID C 60 : 
BIFURCATION APPROACH 

E.E.Tareyeva and T.I.Schelkacheva 

A simple model for the angular dependent interaction between Cqq molecules on the 
face-centered cubic lattice is proposed. The bifurcations of the solutions of the nonlin- 
ear integral equations for orientational distribution functions in the mean field approx- 
imation are analyzed and the orientational phase transition in solid Cqq is described. 
The quantitative results for the orientational phase transition are in accordance with 
experimental data. 

1. As is well known the mean-field approach often brings one to the formulation 
of the broken space symmetry problem in terms of the bifurcation of solutions of 
nonlinear integral equations for distribution functions or order parameters (see, e.g., 
the review [|TJ). In particular, the bifurcation approach was fruitfully used in the 
case of orientational phase transitions in molecular crystals where the intermolecular 
interaction depends not only on the distance between particles but on the molecule 
orientations, too. When the temperature decreases the free rotation of molecules fixed 
on the crystal lattice sites is replaced by rotations about local frozen axes, probably 
different in different sublattices. These orientational trasitions can be conveniently 
described by changes in the n-particle orientational distribution functions f n (u>i, u> n ). 
This kind of problems was treated by use of bifurcation approach, for example, in the 
papers ||, ||, [|, ^j, ||. Here we will follow our papers 0, ||, ||. Let us describe briefly 
the main steps of the method. 

In the mean-field approximation from the first equation of the Bogoliubov hierar- 
chy ]7) for the equilibrium orientational distribution functions or by minimizing the 
orientational free energy one can obtain the following nonlinear integral equation [[J: 

9M) + ^ £ G i J du&ijfai ^i)e 9iK) = 0; (1) 

here Ui are the angles describing the orientation of the molecule, e.g. Euler angles, 
gi{oJi) = ln[^gr^-], fi(uJi) - is the one-particle orientational distribution function for the 
molecule on the site i , the constants Gi are the normalization constants. If we consider 
molecules fixed on the rigid lattice sites then 

®ij(uJi,UJj) = ^2$„ jT (uJ R )u v (u]i)u T (uJj), (2) 

where u u (u>) - is a full system of functions corresponding to the symmetry of the prob- 
lem. In the case of multipole-multipole electrostatic interaction, usually responsible 
for orientational ordering in molecular crystals, we can write in the nearest neighbor 
approximation 

I 2 

$ij(Ui, Uj) = J2 Cl^Ui^Uiv^UirtUj), (3) 
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I\Ui v are the components of the molecule multipole momentum, the coefficients C U)T (uJij) 
are determined by the configuration of the nearest environment, R is the nearest neigh- 
bor distance. 

Assuming that the order is different in not more than k sublattices (k < z), z being 
the number of the nearest neighbors, we will obtain from ([]]) the system of k nonlinear 
integral equations for the functions gi(ui),i = 1, k: 

9i (u) + \Y,G j I ddeP'W C u , T (cj i:j )u lu (u)u lT (u') = 0; (4) 

3=1 v,t 

1 T 2 

A = ' 
kT R 2 ^ 1 ' 

These equations differ from the well known Hammerstein equations |§ only by the 
presence of the factor Gi depending on gi. 

In the case of finite domain of integration when the fixed point principle is valid 
there exists detailed theory for such equations (see 0). We will use the standard 
methods (see e.g.]l0[). In the high temperature region the system @ has only trivial 
solution gi{uj) = 0, corresponding to the orientationally disordered phase. At the 
bifurcation points X a new solutions with broken symmetry appear. Real temperatures 
are described by X a > 0. For X = X a + fi) the functions gf can be written as 
series in integer or fractional powers of fi. These powers are defined by the bifurcation 
equation corresponding to the system (||) or, finally, by the fact wether the integrals 
/ duJUyUrU^ are zero or nonzero. 

The bifurcation points are the eigenvalues of the system of linear integral equations 
of the main order: 

hi(uj) + Xj2 G °j / ^'YC VtT {uij)ui v {u)ui T {u')hj{u') = 0. (5) 

5=1 v,t 

Here hi (u) are the main terms in the expansion of ^ (u) in powers of \x. If / du'Qij (u, a/) = 
0, the equations of the main approximation do not depend on the powers of \x entering 
the expansion of gi(u) . It is easy to see that hiiuj) = J2u h\ui v {uj), so that the definite 
system of nonzero coefficients h\ and, consequently, the definite symmetry of the new 
phase corresponds to each value a. 

Following this way we can obtain in principle all possible broken symmetry phases 
and points of transitions. In fact this program was realized in rather simple cases of 
hydrogen, methane and others. 

However, in the treatise of the complicated problem of the orientational order in 
fullerite we simplify the program from the very beginning. Now our aim is to describe 
the real orientational order seen in the experiments, and to obtain (without using of 
fitting parameters) the transition temperature and the quantitative characteristics of 
the low temperature ordered phase. Our choice of the simplified model interaction is 
based on the maximum utilization of the symmetry properties of the problem and on 
the exploit of the numerical calculations for the energy of the pair of CqO molecules 
with different mutual orientations, performed by other authors. 
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2. First of all, let us recall that the molecule Ceo is the result of cutting off the 
vertices of icosahedron. As result of this cutting off the 20 triangles - the icosahedron 
faces - are transformed into 20 hexagons and 12 (cutted) vertices with 5 outgoing edges 
- in 12 pentagons. The edges of the resulting body are of two types: 30 remaining 
parts of the initial icosahedron edges (the a type) and 5 x 12 = 60 new pentagon edges 
(type b). The detailed calculation of the ground state of an isolated 6*60 molecule shows 
(see, e.g. p6| ) that the hexagons are benzol-type rings with double bonds (and the 
excess of charge) placed on the a-type edges, while the deficiency of charge occurs to 
appear in the centers of hexagons and pentagons. 

The orientational ordering in solid CqO has been a subject of extensive current 
experimental investigations (see, e.g., [H6L HTl |12l IT 



were performed, too 
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; some theoretical researches 
At ambient temperature the 



molecules C 6 rotate almost freely with centers on the face-centered-cubic (fee) lattice 
sites (the space group is Fm3m). 




Figure 1: The sublattices in fee C { 
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When the temperature decreases to T$ ~ 260-fT the first order orientational phase 
transition takes place: the sites of the initial fee lattice become divided between four 
simple cubic sublattices (see fig.l) with own preferable molecular orientation in each 
sublattice (the space group is Pa3). Moreover, the experiments |TT], |T3[ have 



shown that in the ordered state the electron-rich regions (the interpentagon double 
bonds) face the electron-deficient regions of the neighboring 6*60 molecule: the centers 
of pentagons or the centers of hexagons. It was shown that the ratio of the number 
of molecules in those two states is about pp : pn = 60 : 40 at the phase transition 
temperature and increases when the temperature decreases. This remaining two-level 
orientational disorder is usually believed to cause the orientational glass transition at 
Tq ?s 90K now confirmed by various experimental technics (see, e.g., ||24|| ). 

These two minima of the intermolecular angle dependent energy were in fact ob- 
tained through numerical calculations and were shown to be much lower than the 
energies of other mutual orientations of the pair of molecules (see, e.g. 
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In those calculations the previously obtained charge distribution for the isolated C 60 
molecule was taken into account. Usually recent calculations use the intermolec- 
ular potential of Sprik et al. [27j: a sum of 6-12 and Coulomb interactions between 
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60 atoms C and 30 double-bond centers D and between each other: 



fceC(i) jt'eC(2) l^-^W/ \Rkk'J 



k^k',k,k'eC,D {^ K kk'/ \n kk ,J J 

+ E E ^{(^"-(DV S f*- (•> 

fcGD(l) fc'GD(2) I ^-^kk' / \^kk'/ ) k,k'eC,D ^k 1 

Here e = 1.293meV, a C c = 3.4A,<7<7d = 3.5A, crap = 3.6A,g£> = — 0.35e,g<7 = —qp,/2. 
Rigorously speaking to write the equation (HD in our case we are interested in the 
angular part of this complicated interaction represented in terms of multipole-multipole 
interaction of point-like multipoles on the sites of rigid fee lattice with coefficients to 
be calculated from (|J). The general form of this angular part is: 

= C l V)T (uij)ui v (ijJi)ui T (u;j) (7) 

with I = 6,10,12,16,18,... due to the icosahedral molecular symmetry I h (see, e.g., 
0)- 

3. However, we will not use the multipole series per se, but will simplify the problem 
and develop a model orientational interaction exploiting the symmetry properties to 



maximal extend. We will follow the main ideas of the papers |TT|, |21| and use the 
restricted number of allowed orientations instead of free continuous rotations. Let us 
note that a kind of more general discrete model was proposed in Ref. [p9[j . Let us 



take into account in the energy (0) only the orientations with pentagons, hexagons or 
double bonds directed towards 12 nearest neighbors in fee lattice. It is important that 
the C*6o molecule is constructed in such a way that if 6 of its 12 pentagons (or 6 of 
its 20 hexagons) face 6 nearest neighbors double bonds (P and H states of Lapinskas 
et al. then 6 of its 30 interpentagon double bonds (type a) face the remaining 6 

nearest neighbors. 

Now the energy matrix elements can take only three values: Jo ~ the energy of the 
general mutual position, Jp - pentagon versus double bond and Jh - hexagon versus 
double bond. These energies in our model can be compared with those calculated in 
TR [T9|, |2"B"|, 23 as functions of the angular displacements of the molecule at (0,0,0). 



Following and putting J — we obtain from the fig.2(b) of the paper 

J p = -300 K and J H = -110 K. 

The energy matrix elements Jp and Jh connect the states of molecules only in 
the allowed orientations. So, only allowed linear combinations of u\ v enter the energy. 
The theoretical curve in [|T^] makes no difference between the number I of harmonics 
and describes the effect of all of them. So, in the framework of our model calculation 
it is possible to build up the allowed functions using only the harmonics with / = 6: 
we need only their transformation properties. We restrict ourselves to I = 6, however 
the coefficients Jp and Jh are not some of C® T given in (|7|) but effectively take into 
account higher order terms. 
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Let us construct the functions Pi{uS) and Hi(u) explicitly in terms of cubic harmon- 
ics K m = K 6tm ,m = 1,2, ...,13. The states Pi(Hi) have 6 pentagons (hexagons) and 
6 double bonds directed towards 12 nearest neighbors along different [100] axes. All 
functions Pj and Tpv are the sums of K m , invariant under the icosahedral symmetry of 
the molecule (i.e. belonging to the A lg representation of the icosahedral group 7^) if 
icosahedrons are naturally oriented in one of 8 properly chosen coordinate systems. If 
we choose the z-axis along the fifthfold icosahedral axis, the the function A\ (u) can be 

[0: 

A x (u) = a [Y 6fi (uj) + a 5 (r 6i5 (o>) - K 6 _ 5 (a;))], 



written as 



with a = — vTT/5, a 5 = — \/7/\/IT, Yi m {uS) - spherical harmonics (taken as in [ftH]). 
After rotating the icosahedron clockwise about y-axis by the angle = 58.28253 it 
occurs in the position with its three twofold axes parallel to Cartesian axes. This is the 
so called standard P-orientation (following f28j). P\{oj) describes the molecule rotated 
from the standard orientation B about [111] axis through the angle ip = 97.76125° 
(clockwise). After the rotation three fifthfold icosahedron axes occur to be directed 
(almost) precisely along [=pl, 0, =pl], [0, =pl, =Fl], [=f1, =f1, 0], that is in the directions of 
6 nearest neighbors. The functions P 2 (u),P 3 (u) and PAu) are obtained from Pi(uj) 
by subsequent counter-clockwise rotations of the molecule by 90° around z axes. In 
analogous way one can obtain the functions 77j(a;), with the only difference that now 
ip' = 37.76125°. 

If written in the "standard" coordinate frame with Cartesian axes along the cube 
sides these functions have the following explicit form: 



Pi(w) 

p 2 M 

Ps(w) 

PM 

H x {u) 
HAu) 



to 



apKi 
apKi 
apKi 
otpK\ 

a H K x (uj 
a H K x (uj 



ve tne tollowmg explicit torm: 

+(3p[K 8 (u) + K 9 (uj) + K 10 (lo)} + 7p[7iuM + K 12 (u) + K 13 (u)\, 
+f3 P [~K s (uj) + K (u)-K 10 (u>)] + 7p[-^h(^) + K 12 (u)-K 13 (u;)l 
+(3p[K 8 (u)-K 9 (u)-K w (co))+jp[K n ^)-K 12 (u)-K 13 (u)), 
+(3 p [-K s (lo)-K 9 (uj) + K 10 (lo)} + lP [-K ll {u)-K l2 (u) + K 13 (u)\, 
+(3 h [K s {uj) + K 9 {u) + K 10 (u)] + 1h[Kii{u) + K l2 {u) + K 13 {u)\, 
+(3h[-K 8 (uj) + K g (w)-K 10 (u)] + 1h[-Ku{uj) + K X2 (u)-K l3 (oo)\, 
+(3 H [K 8 (u)~K 9 (u)~K 10 (u)}+ lH [K n (oo)~K u (u)-K n (uj)} 1 
+(3 h [-K 8 (uj)-K 9 (u) + K 10 (uj)} + 1h [-K u (uj)-K 12 (uj) + K 13 (lu)}, 

(8) 



with a P = -0.38866; (3 P = 0.31486; 7 P = -0.42877; a H = 0.46588; f3 H = 0.37740; 
jh — 0.34432. The functions are normalized to unity. We use the notations for cubic 
harmonics from |20| (see Appendix). 

So, we suppose in the spirit of the papers 



11, 211 that instead of free rotations 



only discrete jumps between restricted number of orientations is allowed even in the 
high temperature phase. Let us remind that the number of these allowed rotations is 
a rather large one: the molecule has 20 hexagons, 12 pentagons and 30 double bonds. 

4. Let us write now the equations (|) for our model in the case of four sublattices 
of experimentally determined Pa3 structure (see fig. 1): 
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A f 

h^cu) + — / du/[B(w, u')h 2 {uj') + A(u, u')h 3 {uj') + D(lu, oj')h 4 (uj')} = 0, 

47T J 

A /" 

h 2 (cu) + — / + A(o;,a/)^4(^0 + D(u,u')h 3 (u')] = 0, 

47T J 

A r 

h 3 (u) + — / du/fS^u/M^O + A(u,u/)hi(u/) + D(uj,uj')h 2 (uj')} = 0, 
A /" 

/i 4 (w) + — / du'[B(uj, u')h 3 (u') + A(u, u')h 2 {u') + D(u, u/)h 1 (u')] = 0. (9) 

4:71 J 

Here A = 1/T (Jp and Jh are taken in Kelvin) and A(u,u'), B(u,u'), D(u,lu') are 
the sums of interactions over nearest neighbors of the central molecule in the sublattices 
A, B and D (see fig.l), respectively. For example, the sum in the plain perpendicular 
to the x axis (the sites D in fig.l) can be written explicitly in the form: 



D(u,u/) 






2{[(Pi(«) 


+ 


P 4 (u)\ Jp + {H x {u) + Hi(u)] J H ] 


x[P 2 (a/) 


+ 


P 3 (J) + H 2 {J) + H 3 {J)\ 




+ 


P 3 (uj) + H 2 (oj) + H 3 (oj)} 


x[(Pi(u/) 


+ 




+[(AH 


+ 


P 3 (u))Jp+(H 2 (u)+H 3 (u))J h ] 


x[Pi(u/) 


+ 


P a {uj') + H^uj') + H a {uj')} 




+ 


P±{u) + Htiou) + H 4 (uj)} 


x[(P 2 (u/) 


+ 


Ps(u>'))Jp + (H 2 ((J) + H 3 {u'))Jh\) 



(10) 

The bifurcation points for the solutions of initial system are the eigenvalues of the 
system (|9|) while the nontrivial solutions of @ define the symmetry character of the 
new phase. 

It is easy to see that the solutions of the system ([|) can be written as linear 
combinations 

v 

and the system @ is equivalent to the following system of linear algebraic equations: 

A 13 

K + — Y,i A ™ h l + B ruh T 2 + D TV h\\ = 0, 

4?r r=l 

A 13 

h 2 + t- Y,[Ar U hl + B TV K[ + D TU h T 3 ] = 0, 

47T -, 

A 13 

k 3 + 7~ Y,[AtvK + B rvK + D r,h T 2 ] = 0, 
4?T r=l 

A 13 

K + — Y.[ A TvK + B rv hl + D TU h\] = 0, (11) 

4?r T=l 
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v = 1,8, 9,. ..,13. 

Using the explicit form of the matrices A,B,D it is easy to obtain the only nonzero 
elements: 

Ai t i = #i,i = D 1)X = u, 
^8,8 = B 9i9 = £>io,io = v, 

^■11,11 = -Bl2,12 = -Dl3,13 = %■> 
^8,11 = ^.11,8 = -^9,12 — -^12,9 — -^10,13 = -^13,10 = w ■ 

One can write the elements u,v,z and w in terms of the coefficients ap, /3p, 7p, 
&h, Ph, 1h and energies Jp, Jh and obtain the following values: u = 32 ■ 5.046, v = 
32-94.127, z = 32-7.665, w = —32-37.155. The determinant of the algebraic system (|TT|) 
is factorized in 2 x 2 determinants, so that the eigenvalues A a can be easily obtained. 
Among the values X a there are two positive values. The first one Ai = 47r/u corresponds 
to the solution proportional to K\ and is of no interest now. The second A = A& is the 
positive solution of the equation 

1- ±(v + z) + j^- 2 (vz - w 2 ) = 0, (12) 

namely A& = 0.00364f^~ 1 or 7& = 275K. The corresponding nontrivial eigenf unctions 
have h) = 0, and 



hf 


-hf 


= Q{h\ 


-hi) 


? 0; hf 


-hi 1 


= Q(h 8 2 - 


- hi) + 0; 


hf- 


hf = 


Q(hf - 




¥> 0; hf 


-hf 


= Q(hf 


- hf) ± 0; 


hf 


-hf 


= Q(hl 


-hi) 


+ 0; hf 


-hi 2 


= Q(h 9 3 - 


• hi) + 0; 



where 



and 



47T - v\ h 



X b w 



hf + hf = hl + h 8 3 = hf + hi 1 = h 8 2 + h 8 4 

h 



hf + hf = hf + hf = hf + hf = hf + hf 



hf + hf = h\ + hl = hf + hf = h 9 3 + hl = 0. 

If we add the condition for the functions hi(u) to transform one into another under 
the action of the cubic group rotation elements which leave the fee lattice invariant, 
then only three of the coefficients remain to be independent and the functions hi can 
be written in the following form: 

hi(w) = aPi (u) + bHi (lu) + cK x (u) , 

h 2 (u) = aP^ + bH^ + cK^u), (13) 
h 3 (u) = aP^uj) + bHi(uj) + cK^to), 

h A {u) = aP 2 (u) + bH 2 (u) + cK 1 (uj), (14) 

(15) 
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while 

Now 



hx{cj) 
h 2 {u) 
h 3 (u) 



where 



aotp + ban + c = 0. 



r{K % {u) + K 9 {u) + K 10 (u)) + s{K n {u) + K 12 (cu) + K 13 (u)), 
r(-K s (u) + tf 9 (w) - K 10 (u)) + s(-K u (lu) + K 12 (u) - K 13 (u)), 
r(K 8 (u;) - K 9 (u) - if 10 H) + S (ir n (a;) - K 12 (u) - K 13 {u)), 
r(-K s (uj) - K 9 {uo) + K 10 (u)) + s(-K n (u) - K n {u) + K 13 (^)), 



r = a(3 P + s = a^ P + by H , 
s = Qr. 



(16) 



(17) 



Using the numerical values given above we obtain immediately Q = —0.3707, so 
that b : a. We see now that the bifurcation temperature and the ratio of the number of 
molecules in P and H states occur to coincide with the experimental data 0, O, ITH] 



Pp = —t= 0-608; PH = = 0.392. (19) 

a + b a + b 

5. To obtain the remaining unknown coefficient we use the equations of the second 
order in \x. As usually, we will calculate this coefficient from the condition of solvability 
of these equations, i.e. from the condition on right hand sides to be orthogonal to the 
solutions of the homogeneous integral equations. The second order equations have the 
form (let us remind that gi{uS) 



xi(oj) +— du'[A(u,w' 
An J 

x 2 (lu) + ^ I du'[A(u,u' 
An J 

x 3 {u) + — - / duj'[A(uj : uj' 
An J 

xJ^lo) + J dio'[A(uj,Lo' 
where 



= fj,hi(u) + fJ?Xi(u) + A = Ab(l + (jl) 

x 3 (uj') + B(u, uj')x 2 {uj') + D(u, uj')x&(oj' 
x±(uj') + B(u, uj')xi{uj') + D(uj, u')x 3 (uj' 
X\{uj') + B(u, uj')xi(uj') + D(u, uj')x 2 {uj' 
x 2 {uj') + B(lj, u')x 3 (lu') + D(u, uj')xi{uj' 



Ri, 
R2', 

Ra] 



(20) 



Ri 


= hi 


[«>)- 


8tt J 


[ du r [A(u, 




+ B(lo. 


u'){h 2 ) 2 + D{uj. 




R2 


= h 2 


»- 


x b 

8n J 


f du'[A(u, 




+ B(lu. 


J){h x f + D(u, 


wO(to) 2 ]; 


R 3 


= h 3 


[<*)- 


x b 

8n J 


[ du'[A{u, 




+ B(lo. 


uo'){h 4 ) 2 + D{uj. 




R^ 


= h 4 


»- 


x b 

8nJ 


[ du'[A(u, 


u'){h 2 ) 2 


+ B(u. 


oj')(h 3 f + D(u, 


u'){hi?\. (21) 
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The system of nonlinear inhomogeneous equations fl20|) has nontrivial solutions for 
Xi if the right hand sides (|21|) of the system are orthogonal to the solutions of the 
homogeneous equations. In terms of inhomogeneous linear algebraic equations for the 
coefficients xf in the series for Xi in spherical harmonics the analogous condition reads 
that the rank of the system determinant is equal to the rank of the determinant of the 
system examined above. All 16 equations 

J dwRi(w)hj(u}) = (22) 

are identical and give the following equation for the coefficients s,r: 

s 2 + r 2 + ^{r 3 (Miw + M 3 w) + r 2 s[2(M 3 v + M 2 w) + (Mtw + M 3 z)} + 

rs 2 [2(M 3 w + M 2 z) + (M 2 v + M 4 w)} + s 3 (M 2 w + M 4 z)} = 0. (23) 

Here 



Mt =< 8, 9, 10 >= 0.027; M 2 =< 8, 12, 13 >=< 11, 9, 13 >=< 11, 12, 10 >= 0.093; 

M 3 =< 8, 9, 13 >=< 8, 12, 10 >=< 11, 9, 10 >= -0.012; M 4 =< 11, 12, 13 >= -0.011 
- are the nonzero integrals 

< k,m,n >= J duKk(uj)K m (uj)K n (uj). 

Solving equations QTB]), (|18"D and ( p3| ) simultaneously, we obtain finally a = 
—25.7; b = —16.6; c = —2.27. The main term in gi(u>) near the bifurcation point 
has the form: 

g (uj) = -r[aP(uj) + bH(u)) + cK x (u)\, 
T-T b 
T = —- 

6. With this analytical solution near the bifurcation point it is easy to find numerical 
solution of the initial nonlinear integral equations in wider region of temperatures. The 
solution is of the form ( |17]) with the coefficients r(T) and s(T) depending on the 
temperature and goes in both directions from the bifurcation point (see fig.2). The 
solution has the turning point at T t = 290K. The calculation of the orientational part 
of the free energy of the system 

F = J2 j duJ iJ dvj$ij(oJi,oJj)gi(oj i )g j (ojj) + kT^ J ^i^i(^i) In [Angi (a;*)] 

i<j ' i 

was performed for the obtained solutions and the lowest state was identified. The 
solid curves in fig.2 correspond to the stable solution and the dotted curves - to the 
unstable one. At the temperature T c = 288.4^ the first order phase transition takes 
place with the jump of the order parameter. The disordered phase gi(uji) = is stable 
at higher temperature. The ratio of the coefficients r(T) and s(T) slightly depends 
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T[K] 

Figure 2: The coefficients r and s as functions of temperature. 

on the temperature. The relative fractions of molecules in the states P and H in the 
ordered phase at the transition point are: 

M^c) = ~~r = 0.606; p H (T c ) = —— = 0.394. (24) 
a + b a + o 

These values as well as the value of T c are in verygood agreement with the experimental 



data 11, 13, 15 



7. To conclude, we developed a simple model for angle dependent interaction for 
C 60 molecules in the fee cubic lattice. We used rigorous analytic approach based on 
the Lyapunov-Schmidt theory of bifurcation of solutions of nonlinear integral equa- 
tions to treat this model and found the first order orientational phase transition. The 
results for the transition temperature and the characteristics of the ordered phase are 
in quantitative agreement with the experimental data. 

This work was partially supported by Russian Foundation for Basic Researches 
(Grant No. 98-02-16805). 

The authors would like to thank V.N.Ryzhov and V.A.Davydov for useful discus- 
sions. 

Appendix. 

We use the expressions for cubic harmonics from [ED]]. Let us write them explicitly. 



K\ = + y 6 + z 6 )- -(xV + x 4 z 2 + + y 4 x 2 + z 4 x 2 + z 4 y 2 ) + ^x 2 y 2 z 2 ] 

2v27r 2 4 



vl3 

K 8 = ^ Vl05xy{l - I82 2 + 33z 4 
16v27r 



K u = -^—V231xy[3(l - z 2 ) 2 - 16x 2 y 2 



10 



x 2 + y 2 + z 2 = 1. 

The functions K 9 , K 10 and K 12 , K 13 can be obtained by cyclic permutations of the 
coordinates from K 8 and K u , respectively. 
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